PsL Monthly 1993 December/PSL Monthly Shareware CD-ROM (December 1993).iso/prgmming/dos/pascal/p_mat.exe/PMAT.PAS

Pascal/Delphi Source File | 1993-01-23 | 23.5 KB | 895 lines

{
P-Mat - A Turbo Pascal program for Recursive Matrix Algebra

Version: 1.2
Author: Mark Von Tress, Ph.D.
Date: 01/30/93

Copyright(c) Mark Von Tress 1993


DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY
WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS
ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING
FROM USE OF THIS PROGRAM.
}

Unit pmat;

Interface

Type
     fp = ^double;
     ap = Array[1..1] Of double;
     apptr = ^ap;
     app = Array[1..1] Of apptr;

     vmatrixptr = ^vmatrix;
     vmatrix = Object
        r, c : integer;
        Function m( i, j: integer ): double;
        Function mm( i, j: integer ) : fp;
        constructor MakeMatrix( vr, vc: integer );
        Destructor KillVmatrix;
        Procedure Garbage;
        Procedure Show( strng: String );
        Procedure infomatrix( strng: String );

        private
           v,vcheck: ^app;
           nelems : longint;
           onstack : boolean;
           Procedure purgevectors;
           Procedure allocvectors( rr, cc: integer );
     End;

     mstackptr = ^mstack;
     mstack = Object
        constructor InitMatrixStack;
        Destructor KillMatrixStack;
        Procedure inclevel;
        Procedure declevel;
        Function ReturnMat: vmatrixptr;
        Function DecReturn: vmatrixptr;
        Procedure push( Var a: vmatrixptr );
        Procedure nrerror( strng: String );
        Procedure DumpStack;

        private
          nummats, level : integer;
          next : mstackptr;
          mv : vmatrixptr;
          Procedure cleanstack( a: vmatrixptr );
          Function pop: vmatrixptr;
     End;

Var
    dispatch : mstackptr;

Function matequals( Var lop: vmatrixptr; rop: vmatrixptr ) : vmatrixptr;
Function ident( n: integer ): vmatrixptr;

Function Mult( A, B: vmatrixptr ):vmatrixptr;
Function emult( a, b: vmatrixptr ):vmatrixptr;
Function scmult( a: double; b: vmatrixptr ): vmatrixptr;
Function multsc( b: vmatrixptr; a: double ): vmatrixptr;

Function divis( a, b: vmatrixptr ):vmatrixptr;
Function scdivis( a: double; b: vmatrixptr ): vmatrixptr;
Function divissc( b: vmatrixptr; a: double ): vmatrixptr;

Function add( a, b: vmatrixptr ): vmatrixptr;
Function scadd( a: double; b: vmatrixptr ): vmatrixptr;
Function addsc( b: vmatrixptr; a: double ): vmatrixptr;

Function neg( A: vmatrixptr ): vmatrixptr;
Function sub( A, B: vmatrixptr ):vmatrixptr;
Function scsub( A: double; B: vmatrixptr ):vmatrixptr;
Function subsc( B: vmatrixptr; A: double ):vmatrixptr;

Function tran( a: vmatrixptr ): vmatrixptr;
Function sweep( A: vmatrixptr; k1, k2: integer ) : vmatrixptr;
Function inv( Amat: vmatrixptr ): vmatrixptr;

Function fill( rr, cc: integer; d: double ):vmatrixptr;
Function submat( a: vmatrixptr; lr, r, lc, c: integer ): vmatrixptr;
Function cv( a, b: vmatrixptr ): vmatrixptr;
Function ch( a, b: vmatrixptr ): vmatrixptr;
Function vecdiag( a: vmatrixptr ): vmatrixptr;

Function reada( fid: String ): vmatrixptr;
Procedure writea( fid: String; a: vmatrixptr; titlestr: String );

Procedure setwid( wid: integer );
Procedure setdec( decimals: integer );
Function getwid: integer;
Function getdec: integer;

Implementation

{///////////////// vmatrix functions ///////////////}
{$R-}
Procedure vmatrix.allocvectors;
Var i: integer;
Begin
    if c > 8192 then
       dispatch^.nrerror('A matrix cannot have more than 8192 columns');
    getmem( v, r * sizeof( apptr ) );
    For i := 1 To r Do
        getmem( v^[i], c * sizeof( double ) );
End;

Procedure vmatrix.purgevectors;
Var i: integer;
Begin
    If v = Nil Then exit;
    For i := r Downto 1 Do
        freemem( v^[i], c * sizeof( double ) );
    freemem( v, r * sizeof( apptr ) );
End;
{$R+}

constructor vmatrix.MakeMatrix;
Begin
    r := vr;
    c := vc;
    onstack := false;
    if c > 8192 then
       dispatch^.nrerror('A matrix cannot have more than 8192 columns');
    nelems := longint( longint( r ) * longint( c ) ) * longint( sizeof( double ) );
    If nelems < maxAvail - 256 Then allocvectors( r, c )
    Else dispatch^.nrerror( 'cannot allocate matrix' );
    vcheck := v;
End; { MakeMatrix }

Destructor vmatrix.KillVmatrix;
Begin
    purgevectors;
    vcheck := Nil;
End;

Procedure vmatrix.Garbage;
Var errorint : boolean;
Begin
    errorint := false;
    If vcheck <> v Then errorint := true;
    If v = Nil Then errorint := true;
    If r < 1 Then errorint := true;
    If c < 1 Then errorint := true;
    If errorint Then dispatch^.Nrerror( 'garbage' );
End;

{$R-}

Function vmatrix.m;
Begin
    { access a member on the right side of assignment }
    If ( i < 1 ) Or( j < 1 ) Or( i > r ) Or( j > c ) Then
        dispatch^.nrerror( 'm: index out of range' );
    m := v^[i]^[j];
End;

Function vmatrix.mm;
Begin
    { store a matrix element on the left side of assignment }
    If ( i < 1 ) Or( j < 1 ) Or( i > r ) Or( j > c ) Then
        dispatch^.nrerror( 'mm: index out of range' );
    mm := @v^[i]^[j];
End;

{$R+}

Procedure CopySD( Var source, dest: vmatrixptr );
Var
  i,j : integer;
Begin
    source^.Garbage;
    dest^.Garbage;

    If source = dest Then exit;
    If source^.onstack Then
        With dest^ Do Begin
            purgevectors;
            nelems := source^.nelems;
            v := dispatch^.Pop^.v;
            r := source^.r;
            c := source^.c;
            vcheck := v;
            source^.v := Nil;
            dispose( source, killvmatrix );
            exit;
        End;

    If source^.v = dest^.v Then
        With dest^ Do Begin
            r := source^.r;
            c := source^.c;
            allocvectors( r, c );
            vcheck := v;
            nelems := source^.nelems;
        End;

    If ( dest^.r <> source^.r ) Or( dest^.c <> source^.c ) Then
        With dest^ Do Begin
            purgevectors;
            r := source^.r;
            c := source^.c;
            allocvectors( r, c );
            vcheck := v;
            nelems := source^.nelems;
        End;

    For i := 1 To dest^.r Do
        For j := 1 To dest^.c Do
            dest^.mm( i, j )^ := source^.m( i, j );
End;

Function matequals;
     Begin
         rop^.Garbage;
         lop^.Garbage;
         copySD( rop, lop );
         dispatch^.CleanStack( lop );
         matequals := lop;
     End;


{////////////////////// stack functions /////////}


constructor mstack.InitMatrixStack;
Begin
    nummats := 0;
    level := 1;
    mv := Nil;
    getmem( next, sizeof( mstack ) );
    With next^ Do Begin
        level := 0;
        next := Nil;
        mv := Nil;
        nummats := 0;
    End;
End;

Destructor mstack.KillMatrixStack;
Begin
    freemem( next, sizeof( mstack ) );
End;

Procedure mstack.nrerror;
Begin
    writeln( 'Fatal Error: ', strng );
    writeln( 'Exiting to system' );
    halt;
End;

Procedure mstack.inclevel;
Begin
    level := level + 1;
End;

Procedure mstack.declevel;
Begin
    level := level - 1;
    If level <= 0 Then nrerror( 'declevel: decremented too often' );
End;

Function mstack.ReturnMat;
Begin
    ReturnMat := next^.mv;
End;

Function mstack.DecReturn;
Begin
    declevel;
    DecReturn := next^.mv;
End;

Procedure mstack.push;
   Var
      t : mstackptr;
Begin
    a^.Garbage;
    getmem( t, sizeof( mstack ) );
    t^.level := dispatch^.level;
    t^.next := dispatch^.next;
    t^.mv := a;
    a^.onstack := true;
    dispatch^.next := t;
    dispatch^.nummats := dispatch^.nummats + 1;
    t^.nummats := dispatch^.nummats;
End;

Function mstack.pop;
Var t : mstackptr;
    vv : vmatrixptr;
Begin
    If dispatch^.level < 0 Then
        nrerror( 'pop: missmatched inclevel-declevel statments, level < 0' );
    If ( dispatch^.next^.next = Nil ) Or( Nil = dispatch^.next^.mv ) Then
        nrerror( 'pop: popping off the tail' );

    t := dispatch^.next;
    dispatch^.next := dispatch^.next^.next;
    dispatch^.nummats := dispatch^.nummats - 1;
    t^.mv^.Garbage;
    vv := t^.mv;
    freemem( t, sizeof( mstack ) );
    t := Nil;
    vv^.onstack := false;
    pop := vv;
End;

Procedure mstack.cleanstack;
Var vv : vmatrixptr;
Begin
    While ( dispatch^.next^.next <> Nil ) And
            ( dispatch^.next^.level >= dispatch^.level ) Do Begin
        vv := pop;
        If vv <> a Then
            dispose( vv, KillVmatrix );
    End;
End;

Procedure mstack.dumpstack;
Var mm : mstackptr;
Begin
    mm := dispatch;
    writeln( '------------ Stack Dump --------------' );
    While mm <> Nil Do Begin
        With mm^ Do Begin
            writeln( '----------------------------------------------' );
            writeln( 'nummats, level, next, vector: ', nummats, ' ', level, ' ',
                longint( next ), ' ', longint( mv ) );
            If mv <> Nil Then mv^.infomatrix( 'stack matrix information' )
            Else
                writeln;
        End;
        mm := mm^.next;
    End;
End;
{////////////////// matrix functions ///////////////////}

Procedure vmatrix.Show;
   Var i,j,w,d: integer;
   Begin
       w := getwid;
       d := getdec;
       writeln;
       writeln( strng );
       writeln;
       For i := 1 To r Do Begin
           For j := 1 To c Do
               write( m( i, j ): w: d, ' ' );
           writeln;
       End;
       writeln;
   End; { show }

Procedure vmatrix.infomatrix;
   Begin
       writeln;
       writeln( strng );
       writeln;
       writeln( 'r,c,size : ', r, ' ', c, ' ', nelems, ' bytes' );
       writeln( 'v, vcheck: ', longint( v ), ' ', longint( vcheck ) );
   End; { show }

Function mult;
    Var i,j,k: integer;
        C: vmatrixptr;
        s: extended;
    Begin
        a^.Garbage;
        b^.Garbage;
        If A^.c <> B^.r Then
            dispatch^.nrerror( 'mult: matrices do not conform' )
        Else Begin
            new( C, MakeMatrix( a^.r, b^.c ) );
            For i := 1 To A^.r Do
                For j := 1 To B^.c Do Begin
                    s := 0;
                    For k := 1 To B^.r Do
                        s := s + A^.m( i, k ) * B^.m( k, j );
                    C^.mm( i, j )^ := s;
                End;
        End;
        Dispatch^.Push( C );
        Mult := Dispatch^.ReturnMat;
    End; {mult}

Function Emult{( A, B: vmatrixptr): vmatrixptr};
    Var i,j: integer;
        C: vmatrixptr;
    Begin
        a^.Garbage;
        b^.Garbage;
        If ( A^.c <> B^.c ) Or( A^.r <> B^.r ) Then
            dispatch^.Nrerror( 'Emult: matrices do not conform' )
        Else Begin
            new( c, MakeMatrix( A^.r, A^.c ) );
            For i := 1 To a^.r Do
                For j := 1 To a^.c Do
                    C^.mm( i, j )^ := A^.m( i, j ) * B^.m( i, j );
        End;
        Dispatch^.push( c );
        emult := Dispatch^.returnmat;
    End; {emult}

Function scmult{( a: double; b: vmatrixptr): vmatrixptr};
    Var i,j: integer;
        C: vmatrixptr;
    Begin
        b^.Garbage;
        new( c, MakeMatrix( b^.r, b^.c ) );
        For i := 1 To b^.r Do
            For j := 1 To b^.c Do
                C^.mm( i, j )^ := a * B^.m( i, j );
        Dispatch^.push( c );
        scmult := Dispatch^.returnmat;
    End; {scmult}

Function multsc{( b: vmatrixptr; a: double): vmatrixptr};
    Var i,j: integer;
        C: vmatrixptr;
    Begin
        b^.Garbage;
        new( c, MakeMatrix( b^.r, b^.c ) );
        For i := 1 To b^.r Do
            For j := 1 To b^.c Do
                C^.mm( i, j )^ := a * B^.m( i, j );
        Dispatch^.push( c );
        multsc := Dispatch^.returnmat;
    End; {multsc}

Function divis{( A, B: vmatrixptr): vmatrixptr};
    Var i,j: integer;
        C: vmatrixptr;
        dd: double;
    Begin
        a^.Garbage;
        b^.Garbage;
        If ( A^.c <> B^.c ) Or( A^.r <> B^.r ) Then
            dispatch^.Nrerror( 'Emult: matrices do not conform' )
        Else Begin
            new( c, MakeMatrix( b^.r, b^.c ) );
            For i := 1 To a^.r Do
                For j := 1 To a^.c Do Begin
                    dd := b^.m( i, j );
                    If dd = 0.0 Then
                        dispatch^.nrerror( 'DIVIS: zero divisor' );
                    C^.mm( i, j )^ := A^.m( i, j ) / dd;
                End;
        End;
        Dispatch^.push( c );
        divis := Dispatch^.returnmat;
    End; {divis}

Function scdivis{( A : double; B: vmatrixptr): vmatrixptr};
    Var i,j: integer;
        C: vmatrixptr;
        dd: double;
    Begin
        b^.Garbage;
        new( c, MakeMatrix( b^.r, b^.c ) );
        For i := 1 To b^.r Do
            For j := 1 To b^.c Do Begin
                dd := b^.m( i, j );
                If dd = 0.0 Then
                    dispatch^.nrerror( 'SCDIVIS: zero divisor' );
                C^.mm( i, j )^ := A / dd;
            End;
        Dispatch^.push( c );
        scdivis := Dispatch^.returnmat;
    End; {divis}

Function divissc{( A : double; B: vmatrixptr): vmatrixptr};
    Var i,j: integer;
        C: vmatrixptr;
    Begin
        b^.Garbage;
        If a = 0.0 Then
            dispatch^.nrerror( 'DIVISSC: zero divisor' );
        new( c, MakeMatrix( b^.r, b^.c ) );
        For i := 1 To b^.r Do
            For j := 1 To b^.c Do Begin
                C^.mm( i, j )^ := b^.m( i, j ) / A;
            End;
        Dispatch^.push( c );
        divissc := Dispatch^.returnmat;
    End; {divissc}

Function ADD{( A, B: vmatrixptr): vmatrixptr};
    Var i,j: integer;
        C: vmatrixptr;
Begin
    a^.Garbage;
    b^.Garbage;
    If ( A^.c <> B^.c ) Or( A^.r <> B^.r ) Then
        dispatch^.Nrerror( 'add: matrices do not conform' )
    Else Begin
        new( c, MakeMatrix( A^.r, A^.c ) );
        For i := 1 To a^.r Do
            For j := 1 To a^.c Do
                C^.mm( i, j )^ := A^.m( i, j ) + B^.m( i, j );
    End;
    Dispatch^.push( c );
    ADD := Dispatch^.returnmat;
End; {add}

Function scADD{( A : double; B: vmatrixptr): vmatrixptr};
    Var i,j: integer;
        C: vmatrixptr;
Begin
    b^.Garbage;
    new( c, MakeMatrix( b^.r, b^.c ) );
    For i := 1 To b^.r Do
        For j := 1 To b^.c Do
            C^.mm( i, j )^ := A + B^.m( i, j );
    Dispatch^.push( c );
    scADD := Dispatch^.returnmat;
End; {scadd}

Function ADDSC{(  B : vmatrixptr; A : double): vmatrixptr};
    Var i,j: integer;
        C: vmatrixptr;
Begin
    b^.Garbage;
    new( c, MakeMatrix( b^.r, b^.c ) );
    For i := 1 To b^.r Do
        For j := 1 To b^.c Do
            C^.mm( i, j )^ := A + B^.m( i, j );
    Dispatch^.push( c );
    ADDsc := Dispatch^.returnmat;
End; {addsc}

Function neg;
  Var i,j: integer;
  B: vmatrixptr;
  Begin
      a^.garbage;
      new( b, MakeMatrix( a^.r, a^.c ) );
      For i := 1 To a^.r Do
          For j := 1 To a^.c Do
              B^.mm( i, j )^ := - A^.m( i, j );
      dispatch^.Push( b );
      neg := dispatch^.returnmat;
  End; {negate}
Function sub{( A, B: vmatrixptr): vmatrixptr};
    Var i,j: integer;
        C: vmatrixptr;
Begin
    a^.Garbage;
    b^.Garbage;
    If ( A^.c <> B^.c ) Or( A^.r <> B^.r ) Then
        dispatch^.Nrerror( 'sub: matrices do not conform' )
    Else Begin
        new( c, MakeMatrix( A^.r, A^.c ) );
        For i := 1 To a^.r Do
            For j := 1 To a^.c Do
                C^.mm( i, j )^ := A^.m( i, j ) - B^.m( i, j );
    End;
    Dispatch^.push( c );
    sub := Dispatch^.returnmat;
End; {sub}

Function scsub{ (a: double; b:vmatrixptr): vmatrixptr;};
    Var i,j: integer;
        C: vmatrixptr;
    Begin
        b^.Garbage;
        new( c, MakeMatrix( b^.r, b^.c ) );
        For i := 1 To b^.r Do
            For j := 1 To b^.c Do
                C^.mm( i, j )^ := A - B^.m( i, j );
        dispatch^.Push( c );
        scSub := dispatch^.returnmat;
    End; { scsub }

Function subsc{(b: vmatrixptr; a : double): vmatrixptr;};
    Var i,j: integer;
        C: vmatrixptr;
    Begin
        b^.Garbage;
        new( c, MakeMatrix( b^.r, b^.c ) );
        For i := 1 To b^.r Do
            For j := 1 To b^.c Do
                C^.mm( i, j )^ := B^.m( i, j ) - A;
        dispatch^.Push( c );
        Subsc := dispatch^.returnmat;
    End; { subsc }

Function tran;
    Var i,j: integer;
        c : vmatrixptr;
Begin
    a^.garbage;
    new( c, MakeMatrix( A^.c, A^.r ) );
    For i := 1 To A^.r Do
        For j := 1 To A^.c Do
            c^.mm( j, i )^ := A^.m( i, j );
    dispatch^.Push( c );
    tran := dispatch^.returnmat;
End; {transpose}


Function sweep;
{
     SUBROUTINE TO SWEEP THE NXN MATRIX A ON ITS K1 THRU K2
     DIACONAL ELEMENTS (SWP(K)SWP(K)A=A)

     INPUT : A Matrix to be swept
             K1,K2: elements to sweep; if k1=1 and k2= n then gets inverse
}
Const Eps = 1.0E-8;
Var I,J,k,n,it: integer;
     d,ss: extended;
Begin
    a^.garbage;
    If A^.r <> A^.c Then
        dispatch^.nrerror( 'sweep: Amat not square in Sweep: ' );

    n := A^.r;
    If K2 < K1 Then Begin
        k := k1; k1 := k2; k2 := k;
    End;
    For K := K1 To K2 Do Begin
        If Abs( A^.m( K, K ) ) < eps Then Begin
            For it := 1 To n Do Begin
                a^.mm( it, k )^ := 0;
                a^.mm( k, it )^ := 0;
            End;
        End
        Else Begin
            D := 1.0 / A^.m( K, K );
            A^.mm( K, K )^ := 1.0;
            For I := 1 To N Do A^.mm( K, I )^ := D * A^.m( K, i );
            For J := 1 To N Do
                If J <> K Then A^.mm( J, k )^ :=
                        - D * A^.m( J, K );
            For J := 1 To N Do
                If J <> K Then
                    For I := 1 To N Do Begin
                        ss := A^.m( J, I );
                        If I <> K Then A^.mm( J, I )^ :=
                                ss + (1.0 / D) * A^.m( J, K ) * A^.m( K, I );
                    End;
        End;
    End;
    dispatch^.push( a );
    sweep := dispatch^.returnmat;
End; {sweep}

Function inv;
   Var AMat1: vmatrixptr;
   Begin
       dispatch^.inclevel;
       amat^.garbage;
       If AMat^.r <> AMat^.c Then
           dispatch^.nrerror( 'Matrix not square in Invert ' );
       new( amat1, Makematrix( amat^.r, amat^.r ) );
       amat1 := matequals( amat1, amat );
       amat1^.show( 'amat1' );
       AMat1 := matequals( amat1, Sweep( amat1, 1, amat1^.r ) );
       dispatch^.push( amat1 );
       inv := dispatch^.Decreturn;
   End;  {Inv}

{//////////////////// patterned matrices ///////////////////}

Function ident;
  Var i,j: integer;
      a: vmatrixptr;
Begin
    new( a, MakeMatrix( n, n ) );
    For i := 1 To n Do
        For j := 1 To n Do
            If i = j Then a^.mm( i, i )^ := 1.0
            Else a^.mm( i, j )^ := 0;
    dispatch^.Push( a );
    ident := dispatch^.ReturnMat;
End; { ident }


Function fill;
  Var i,j: integer;
      c : vmatrixptr;
  Begin
      If ( rr < 1 ) Or( cc < 1 ) Then
          dispatch^.nrerror( 'fill: negative rows or columns' );
      new( c, makematrix( rr, cc ) )    ;
      For i := 1 To rr Do
          For j := 1 To cc Do
              c^.mm( i, j )^ := d;
      dispatch^.push( c );
      fill := dispatch^.returnmat;
  End;

Function submat;
   Var i,j,rc,cc: integer;
       b : vmatrixptr;
Begin
    a^.Garbage;
    If ( lr < 1 ) Or( r > a^.r ) Or( lc < 1 ) Or( c > a^.c ) Then dispatch^.nrerror( 'SUBMAT: index out of range' );
    If ( lr > r ) Or( lc > c ) Then
        dispatch^.nrerror( 'SUBMAT: indexes out of order, lc>c or lr>r' );

    new( b, makematrix( r - lr + 1, c - lc + 1 ) );

    rc := b^.r;
    cc := b^.c;
    For i := 1 To rc Do
        For j := 1 To cc Do
            b^.mm( i, j )^ := a^.m( lr - 1 + i, lc - 1 + j );

    dispatch^.push( b );
    submat := dispatch^.ReturnMat;
End;

Function ch;
  Var i,j,k: integer;
      c : vmatrixptr;
Begin
    { horizontal concatination }
    a^.Garbage;
    b^.Garbage;
    If a^.r <> b^.r Then
        dispatch^.nrerror( 'CH: matrices have different number of rows' );

    k := a^.c + b^.c;
    new( c, makematrix( a^.r, k ) );

    For i := 1 To a^.r Do
        For j := 1 To k Do
            If j <= a^.c Then c^.mm( i, j )^ := a^.m( i, j )
            Else c^.mm( i, j )^ := b^.m( i, j - a^.c );

    dispatch^.push( c );
    ch := dispatch^.ReturnMat;
End;     { ch }

Function cv;
  Var i,j,k : integer;
          c : vmatrixptr;
Begin
    { vertical concatination }
    a^.Garbage;
    b^.Garbage;
    If a^.c <> b^.c Then
        dispatch^.nrerror( 'CV: matrices have different number of rows' );

    k := a^.r + b^.r;
    new( c, MakeMatrix( k, a^.c ) );

    For i := 1 To a^.c Do
        For j := 1 To k Do
            If j <= a^.r Then c^.mm( j, i )^ := a^.m( j, i )
            Else c^.mm( j, i )^ := b^.m( j - a^.r, i );

    dispatch^.push( c );
    cv := dispatch^.ReturnMat;
End;     { cv }

Function vecdiag;
Var  i,j,r: integer;
     b : vmatrixptr;
     t : double;
Begin
    { insert a vector into the diagonal of I }
    a^.Garbage;
    r := 0;
    If a^.r = 1 Then r := a^.c;
    If a^.c = 1 Then r := a^.r;
    If r = 0 Then dispatch^.nrerror( 'vecdiag: a is not a vector' );
    new( b, MakeMatrix( r, r ) );

    For i := 1 To r Do Begin
        If a^.r = 1 Then t := a^.m( 1, i )
        Else  t := a^.m( i, 1 );
        For j := 1 To r Do
            If i = j Then b^.mm( i, j )^ := t
            Else b^.mm( i, j )^ := 0.0
    End;

    dispatch^.push( b );
    vecdiag := dispatch^.ReturnMat;
End;     {vecdiag}

Function FileExists( FileName: String )
                                : Boolean;
{ Returns True if file exists; otherwise,
  it returns False. Closes the file if
  it exists. }
Var
  f: File;
Begin
    {$I-}
    Assign( f, FileName );
    Reset( f );
    Close( f );
    {$I+}
    FileExists := (IOResult = 0) And
        (FileName <> '');
End;  { FileExists }

Function ReadA;
Var i,j  : integer;
    b : vmatrixptr;
    f : text;
    titlestr: String;
Begin
    If Not FileExists( fid ) Then Begin
        writeln( 'trying to open : ', fid );
        dispatch^.nrerror( 'reada: file does not exist' );
    End;

    assign( f, fid );
    reset( f );

    writeln( 'reading: ', fid );
    readln( f, titlestr );
    writeln( 'title:' );
    writeln( titlestr );
    read( f, i, j );
    new( b, makematrix( i, j ) );
    For i := 1 To b^.r Do
        For j := 1 To b^.c Do
            read( f, b^.mm( i, j )^ );
    close( f );
    dispatch^.push( b );
    reada := dispatch^.returnmat;
End; { reada }

Procedure writea;
Var i,j,w,d  : integer;
    f :  text;
Begin
    assign( f, fid );
    rewrite( f );

    w := getwid;
    d := getdec;
    writeln( f, titlestr );
    writeln( f, a^.r, ' ', a^.c );
    For i := 1 To a^.r Do Begin
        For j := 1 To a^.c Do
            write( f, a^.m( i, j ): w: d, ' ' );
        writeln( f );
    End;
    close( f );
End; { reada }

Var
   printdec, printwid : integer;



Function getwid;
Begin
    getwid := printwid;
End;

Function getdec;
Begin
    getdec := printdec;
End;

Procedure setwid;
Begin
    If wid > 0 Then printwid := wid;
End;

Procedure  setdec;
Begin
    If ( decimals > 0 ) And( decimals <= printwid ) Then
        printdec := decimals;
End;

Begin
    setwid( 6 );
    setdec( 2 );
    new( dispatch, InitMatrixStack );
End.